/**
 * 
 */
package ua.lviv.franko.solvers.stressanalysis.twoDimention.combined;

import java.util.ArrayList;
import java.util.Collections;
import java.util.HashSet;
import java.util.List;
import java.util.Set;

import org.jfree.data.xy.XYSeries;
import org.jfree.data.xy.XYSeriesCollection;
import org.jzy3d.maths.Range;
import org.jzy3d.plot3d.builder.Mapper;

import ua.lviv.franko.integration.ConstantFunction1D;
import ua.lviv.franko.integration.ConstantFunction2D;
import ua.lviv.franko.integration.FunctionWithoutDerivatives1D;
import ua.lviv.franko.integration.FunctionWithoutDerivatives2D;
import ua.lviv.franko.integration.Gauss2D;
import ua.lviv.franko.integration.IFunction1D;
import ua.lviv.franko.integration.IFunction2D;
import ua.lviv.franko.matrix.hash.lu.LUFactorization;
import ua.lviv.franko.solvers.stressanalysis.oneDimention.AbstractElement1D;
import ua.lviv.franko.solvers.stressanalysis.oneDimention.EquationParameters;
import ua.lviv.franko.solvers.stressanalysis.oneDimention.FEElementLinear;
import ua.lviv.franko.solvers.stressanalysis.oneDimention.TimoshenkoMembraneSolver;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.Solver;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.AbstractStressElement;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.fem.LinearRectangularElement;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.savestructures.bounderyconditions.BounderyConditionLinear;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.savestructures.bounderyconditions.BounderyConditionStressInterface;
import ua.lviv.franko.solvers.stressanalysis.twoDimention.savestructures.bounderyconditions.BounderyConditionType;
import ua.lviv.franko.triangulation.SaveStructures.Node;
import ua.lviv.franko.triangulation.SaveStructures.NodePair;
import ua.lviv.franko.triangulation.SaveStructures.Rectangle;
import ua.lviv.franko.visualisation.NewOkCancelDialog;
import ua.lviv.franko.visualisation.Visualise3D;

/**
 * @author marko
 * 
 */
public class DomainDecompositionDirihletNauman {
	public static final double					EPS			= 0.000000000000001;
	protected TimoshenkoMembraneSolver			membraneSolver;
	protected Solver							solver;
	protected double							epsPerc;

	protected BounderyConditionStressInterface	bounderyConditionsTop;
	protected double							yTopValue	= 1.0;

	public DomainDecompositionDirihletNauman(TimoshenkoMembraneSolver membraneSolver, Solver solver, double epsPerc, BounderyConditionStressInterface bounderyConditionsTop) {
		super();
		this.membraneSolver = membraneSolver;
		this.solver = solver;
		this.epsPerc = epsPerc;
		this.bounderyConditionsTop = bounderyConditionsTop;
		solver.generate();
		solver.applyBounderyConditions(solver.getBounderyConditions());
		membraneSolver.generate();
		membraneSolver.setBothDirihlet(new double[] { 0.0, 0.0, 0.0 }, new double[] { 0.0, 0.0, 0.0 });
	}

	public void solve() {
		final double[] results2D = new double[solver.getB().length];
		final double[] membraneVector1D = new double[membraneSolver.getB().length];
		final EquationParameters old = membraneSolver.getElementsList().get(0).getEquationParameters();
		EquationParameters eParameters = new EquationParameters();
		eParameters.setA1(old.getA1());
		eParameters.setE(old.getE());
		eParameters.setF1(old.getF1());
		eParameters.setF3(old.getF3());
		eParameters.setG(old.getG());
		eParameters.setH(old.getH());
		eParameters.setK1(old.getK1());
		eParameters.setkPrime(old.getkPrime());
		eParameters.setSigma11Plus(old.getSigma11Plus());
		eParameters.setSigma13Plus(old.getSigma13Plus());
		eParameters.setSigma33Plus(old.getSigma33Plus());
		eParameters.setSigma11Minus(old.getSigma11Minus());
		eParameters.setSigma13Minus(old.getSigma13Minus());
		eParameters.setSigma33Minus(old.getSigma33Minus());
		final List<AbstractStressElement> stressElements = solver.getElementsList();

		LUFactorization membraneLU = new LUFactorization(membraneSolver.getMainMatrix());

		// LUFactorization stressLU = new LUFactorization(solver.getMainMatrix());
		int iterations = 0;
		while (true) {
			iterations++;
			//2d
			IFunction2D topStressConditions1 = new FunctionWithoutDerivatives2D() {
				// u1 = -w
				@Override
				public double calculate(double x, double y) {
					List<AbstractElement1D> element1ds = membraneSolver.getElementsList();
					double res = 0.0;
					for (int i = 0; i < element1ds.size(); i++) {
						double[] tmp = element1ds.get(i).calculate(x, membraneVector1D);
						res += tmp[1];
						if (Math.abs(res)>EPS) {
							break;
						}
					}
					return res;
				}
			};

			IFunction2D topStressConditions2 = new FunctionWithoutDerivatives2D() {
				// u2 = -v1 + gama*h/2
				@Override
				public double calculate(double x, double y) {
					List<AbstractElement1D> element1ds = membraneSolver.getElementsList();
					double hDivTwo = old.getH() / 2.0;
					double res = 0.0;
					for (int i = 0; i < element1ds.size(); i++) {
						double[] tmp = element1ds.get(i).calculate(x, membraneVector1D);
						res += tmp[0] + hDivTwo * tmp[2];
						if (Math.abs(res)>EPS) {
							break;
						}
					}
					return res;
				}
			};
			bounderyConditionsTop.setBounderyValue1(topStressConditions1);
			bounderyConditionsTop.setBounderyValue2(topStressConditions2);
			List<BounderyConditionStressInterface> bndrs = new ArrayList<BounderyConditionStressInterface>();
			bndrs.add(bounderyConditionsTop);
			solver.applyBounderyConditions(bndrs);
			double[] mainProblemResults = solver.solve();
			for (int i = 0; i < mainProblemResults.length; i++) {
				results2D[i] = mainProblemResults[i];
			}
			
			//end 2d
			eParameters.setSigma11Minus(new FunctionWithoutDerivatives1D() {
				@Override
				public double calculate(double x) {
					double val = 0.0;
					int i = 0;
					while (val == 0.0 && i < stressElements.size()) {
						val += -stressElements.get(i).calculateSigma22(x, yTopValue, results2D);
						i++;
					}
					return val;
				}
			});
			eParameters.setSigma13Minus(new FunctionWithoutDerivatives1D() {
				@Override
				public double calculate(double x) {
					double val = 0.0;
					int i = 0;
					while (val == 0.0 && i < stressElements.size()) {
						val += -stressElements.get(i).calculateSigma12(x, yTopValue, results2D);
						i++;
					}
					return val;
				}
			});
			double[] membraneVector = membraneSolver.getUpdatedBVector(eParameters);
			TimoshenkoMembraneSolver.setBothDirihlet(membraneSolver.getMainMatrix(), membraneVector, new double[] { 0.0, 0.0, 0.0 }, new double[] { 0.0, 0.0, 0.0 });
			final double[] membraneResults = membraneLU.calculate(membraneVector);
			for (int i = 0; i < membraneResults.length; i++) {
				membraneVector1D[i] = membraneResults[i];
			}
			XYSeries u1 = new XYSeries("result1", true, true);
			XYSeries u2 = new XYSeries("result2", true, true);

			int NNN = 20;
			double a = 0;
			double b = 1;
			double hStep = membraneSolver.getElementsList().get(0).getRight() - membraneSolver.getElementsList().get(0).getLeft();
			AbstractElement1D el = null;
			for (int j = 0; j < membraneSolver.getElementsList().size(); j++) {
				el = membraneSolver.getElementsList().get(j);
				double hhh = hStep / NNN;
				for (int k = 0; k < NNN; k++) {
					double[] tmp = el.calculate(a + hStep * j + k * hhh, membraneResults);
					double x = a + hStep * j + k * hhh;
					u1.add(x, topStressConditions1.calculate(x, 0));
					u2.add(x, topStressConditions2.calculate(x, 0));
//					u1.add(a + hStep * j + k * hhh, tmp[1]);
//					u2.add(a + hStep * j + k * hhh, tmp[0] + old.getH() / 2.0 * tmp[2]);
				}
			}
			double[] tmp = el.calculate(b, membraneResults);
			u1.add(b, tmp[1]);
			u2.add(b, tmp[0] - old.getH() / 2.0 * tmp[2]);
			
			XYSeries sigma1 = new XYSeries("sigma1", true, true);
			XYSeries sigma2 = new XYSeries("sigma2", true, true);
			
			double myH = (b-a)/100.0;
			for (int i = 0; i <= 100; i++) {
				sigma1.add(a+myH*i, eParameters.getSigma11Minus().calculate(a+myH*i));
				sigma2.add(a+myH*i, eParameters.getSigma13Minus().calculate(a+myH*i));
			}
			

			XYSeriesCollection col = new XYSeriesCollection();
//			col.addSeries(sigma1);
//			col.addSeries(sigma2);
			col.addSeries(u1);
			col.addSeries(u2);
			NewOkCancelDialog d = new NewOkCancelDialog(null, false, col);
			d.setLocationRelativeTo(null);
			d.setVisible(true);
			
			
//			Visualise3D.show(new Mapper() {
//				@Override
//				public double f(double x, double y) {
//					double val = 0.0;
//					int i = 0;
//					while (val == 0.0 && i < stressElements.size()) {
//						val += stressElements.get(i).calculateU(x, y, results2D);
//						i++;
//					}
//					return val;
//				}
//			}, new Range(0, 1), new Range(0, 1));
			

			if (iterations == 10) {

				Visualise3D.show(new Mapper() {
					@Override
					public double f(double x, double y) {
						double val = 0.0;
						int i = 0;
						while (val == 0.0 && i < stressElements.size()) {
							val += stressElements.get(i).calculateU(x, y, results2D);
							if (Math.abs(val)>EPS) {
								break;
							}
							i++;
						}
						return val;
					}
				}, new Range(0, 1), new Range(0, 1));
				
				Visualise3D.show(new Mapper() {
					@Override
					public double f(double x, double y) {
						double val = 0.0;
						int i = 0;
						while (val == 0.0 && i < stressElements.size()) {
							val += stressElements.get(i).calculateV(x, y, results2D);
							if (Math.abs(val)>EPS) {
								break;
							}
							i++;
						}
						return val;
					}
				}, new Range(0, 1), new Range(0, 1));

				// XYSeries s1 = new XYSeries("result1", true, true);
				// XYSeries s2 = new XYSeries("result2", true, true);
				// XYSeries s3 = new XYSeries("result3", true, true);

				// int NNN = 20;
				// double a = 0;
				// double b = 1;
				// double hStep = membraneSolver.getElementsList().get(0).getRight() - membraneSolver.getElementsList().get(0).getLeft();
				// AbstractElement1D el = null;
				// for (int j = 0; j < membraneSolver.getElementsList().size(); j++) {
				// el = membraneSolver.getElementsList().get(j);
				// double hhh = hStep / NNN;
				// for (int k = 0; k < NNN; k++) {
				// double[] tmp = el.calculate(a + hStep * j + k * hhh, membraneResults);
				// s1.add(a + hStep * j + k * hhh, tmp[0]);
				// s2.add(a + hStep * j + k * hhh, tmp[1]);
				// s3.add(a + hStep * j + k * hhh, tmp[2]);
				// }
				// }
				// double[] tmp = el.calculate(b, membraneResults);
				// s1.add(b, tmp[0]);
				// s2.add(b, tmp[1]);
				// s3.add(b, tmp[2]);
				//
				// XYSeriesCollection col = new XYSeriesCollection();
				// col.addSeries(s1);
				// col.addSeries(s2);
				// col.addSeries(s3);
				// NewOkCancelDialog d = new NewOkCancelDialog(null, false, col);
				// d.setLocationRelativeTo(null);
				// d.setVisible(true);

				break;
			}

		}
	}

	/**
	 * @param args
	 */
	public static void main(String[] args) {
		IFunction2D f1 = new ConstantFunction2D(1.0);
		IFunction2D f3 = new ConstantFunction2D(1.0);
		IFunction1D k1 = new ConstantFunction1D(0);
		IFunction1D a1 = new ConstantFunction1D(1);
		IFunction1D sigma11Plus = new ConstantFunction1D(0);
		IFunction1D sigma13Minus = new ConstantFunction1D(0);
		IFunction1D sigma11Minus = new ConstantFunction1D(0);
		IFunction1D sigma33Minus = new ConstantFunction1D(0);
		IFunction1D sigma33Plus = new ConstantFunction1D(0);
		IFunction1D sigma13Plus = new ConstantFunction1D(0);
		double v = 0.37;
		double e = 80000;
		double v2 = v;
		double e2 = e/4.0;
		double h = 0.5;

		EquationParameters eqP = new EquationParameters(a1, k1, e, v, h, f1, f3, sigma13Plus, sigma13Minus, sigma11Plus, sigma11Minus, sigma33Plus, sigma33Minus);
		double a = 0.0;
		double b = 1.0;
		int n = 15;
		double hStep = (b - a) / (double) n;
		ArrayList<AbstractElement1D> elements = new ArrayList<AbstractElement1D>();
		for (int i = 0; i < n; i++) {
			elements.add(new FEElementLinear(eqP, i, a + i * hStep, a + (i + 1) * hStep));
			// elements.add(new BubbleElement(eqP, i, a + i * hStep, a + (i + 1) * hStep));
			// elements.add(new EnrichedLinearElement(eqP, i, a + i * hStep, a + (i+1) * hStep, 5, 1, new Gaussian(), 1.5));
		}

		TimoshenkoMembraneSolver s = new TimoshenkoMembraneSolver(elements);

		double ax = -0.0, bx = 1.0, ay = -0.0, by = 1.0;

		int n2 = 10, m = 10;

		double hx = (bx - ax) / (double) n2, hy = (by - ay) / (double) m;

		// double k12d = e / (1 - v * v);
		// double k22d = e * v / (1 - v * v);
		// double k32d = e / (1 - v * v) * ((1 - v) / 2.0);

		final List<AbstractStressElement> elems = new ArrayList<AbstractStressElement>();
		Set<Node> tmpNodesSet = new HashSet<Node>();

		for (int i = 0; i < m; i++) {
			for (int j = 0; j < n2; j++) {
				Rectangle r = new Rectangle();
				r.A = new Node(ax + j * hx, ay, false, i * (n2 + 1) + j);
				r.B = new Node(ax + (j + 1) * hx, ay, false, i * (n2 + 1) + j + 1);
				r.C = new Node(ax + (j + 1) * hx, ay + hy, false, (i + 1) * (n2 + 1) + j + 1);
				r.D = new Node(ax + j * hx, ay + hy, false, (i + 1) * (n2 + 1) + j);
				tmpNodesSet.add(r.A);
				tmpNodesSet.add(r.B);
				tmpNodesSet.add(r.C);
				tmpNodesSet.add(r.D);

				// LinearRectangularElement elem = new LinearRectangularElement(r.A, r.B, r.C, r.D, k12d, k22d, k32d, f1, f3);
				LinearRectangularElement elem = new LinearRectangularElement(r.A, r.B, r.C, r.D, e2, v2, f1, f3);
				elems.add(elem);
			}
			ay += hy;
		}
		ay = -0.0;
		List<Node> nodes = new ArrayList<Node>();
		nodes.addAll(tmpNodesSet);
		Collections.sort(nodes);
		List<BounderyConditionStressInterface> bounderies = new ArrayList<BounderyConditionStressInterface>();
		List<NodePair> leftNodePairs = new ArrayList<NodePair>();
		List<NodePair> rightNodePairs = new ArrayList<NodePair>();
		List<NodePair> topNodePairs = new ArrayList<NodePair>();
		List<NodePair> bottomNodePairs = new ArrayList<NodePair>();
		for (AbstractStressElement element : elems) {
			List<Integer> positions = element.getPos();
			Node A = nodes.get(positions.get(0));
			Node B = nodes.get(positions.get(1));
			Node C = nodes.get(positions.get(2));
			Node D = nodes.get(positions.get(3));
			// left
			List<Node> tmpNodes = new ArrayList<Node>();
			if (Math.abs(A.x - ax) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.x - ax) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.x - ax) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.x - ax) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				leftNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}

			// right
			tmpNodes.clear();
			if (Math.abs(A.x - bx) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.x - bx) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.x - bx) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.x - bx) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				rightNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}

			// top
			tmpNodes.clear();
			if (Math.abs(A.y - by) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.y - by) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.y - by) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.y - by) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				topNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}

			// bottom
			tmpNodes.clear();
			if (Math.abs(A.y - ay) < EPS) {
				tmpNodes.add(A);
			}
			if (Math.abs(B.y - ay) < EPS) {
				tmpNodes.add(B);
			}
			if (Math.abs(C.y - ay) < EPS) {
				tmpNodes.add(C);
			}
			if (Math.abs(D.y - ay) < EPS) {
				tmpNodes.add(D);
			}

			if (tmpNodes.size() == 2) {
				bottomNodePairs.add(new NodePair(tmpNodes.get(0), tmpNodes.get(1)));
			}
		}
		Gauss2D gauss2d = new Gauss2D(7);

		bounderies.add(new BounderyConditionLinear(BounderyConditionType.DIRIHLET_TYPE, new ConstantFunction2D(0.0), new ConstantFunction2D(0.0), leftNodePairs, gauss2d));

		bounderies.add(new BounderyConditionLinear(BounderyConditionType.DIRIHLET_TYPE, new ConstantFunction2D(0.0), new ConstantFunction2D(0.0), rightNodePairs, gauss2d));

		Solver solver = new Solver((n2 + 1) * (m + 1), elems, bounderies);
		List<AbstractStressElement> nearBounderyElements = new ArrayList<AbstractStressElement>();
		for (AbstractStressElement abstractStressElement : elems) {
			List<Integer> nodeNumbers = abstractStressElement.getNodesNumbers();
			for (Integer index : nodeNumbers) {
				if (Math.abs(nodes.get(index).y - by) < EPS) {
					nearBounderyElements.add(abstractStressElement);
					break;
				}
			}
		}

		DomainDecompositionDirihletNauman decomposition = new DomainDecompositionDirihletNauman(s, solver, 0.001, new BounderyConditionLinear(BounderyConditionType.DIRIHLET_TYPE,
				new ConstantFunction2D(0.0), new ConstantFunction2D(0.0), topNodePairs, gauss2d));
		decomposition.solve();

	}

}
